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This is a progress report on the research project "Numerical Solutions 
of Navier- Stokes Equations for a Butler Wing." The period of performance on 
this research was January 1 through August 31, 1985. The work is supported 
by the NASA/Langley Research Center through Cooperative Agreement NCC1-68, 
and monitored by Dr. Robert E. Smith, Jr., of the Analysis and Computation 
Division (Computer Applications Branch), NASA/Langley, MS/125. 


ABSTRACT 


The flow field is simulated on the surface of a. given delta wing 
(Butler wing) at zero incident in a uniform stream. The simulation is done 
by integrating a set of flow field equations. This set of equations governs 
the unsteady, viscous, compressible, heat conducting flow of an ideal gas. 
The equations are written in curvilinear coordinates so that the wing 
surface is represented accurately. These equations are solved by the finite 
difference method, and results obtained for high-speed freestream conditions 
are compared with theoretical and experimental results. 

In the present study, the Kavier-Stokes equations are solved numerical- 
ly. These equations are unsteady, compressible, viscous, and three-dimen- 
.sional without neglecting any terms. The time dependency of the governing 
equations allows the solution to progress naturally for an arbitrary initial 
guess to an asymptotic steady state, if one exists. The equations are 
transformed from physical coordinates to the computational coordinates, 
allowing the solution of the governing equations in a rectangular parallele- 
piped domain. The equations are solved by the i<1acCorimac.k time-split tech- 
nique which is vectorized and programmed to run on the CDC VPS 32 computer. 
The codes are written in 32-bit (half word) FORTRAN, which provide: an ap- 
proximate factor of two decrease in computational time and doubles the 
memory size compared to the 64-bit word size. 
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NUMERICAL SOLUTIONS OF NAV IE R- STOKES EQUATIONS 
FOR A BUTLER WING 

By 

J. S. Abolhassani 1 and S. N. Tiwari 2 
1. INTRODUCTION 

The Butler wing is a delta wing which was proposed by D. S. Butler 
(Ref. 4). The planform of the body is an isosceles triangle, and the lead- 
ing edges of the wing are laid along the Mach lines of the undisturbed 
stream. For the first 20% of the length, the body is a right circular cone. 
The remainder of the body has elliptic sections which become more eccentric 
as the sharp trailing edges are approached. 

0. S. Butler has compared experimental results for surface pressure 
with theoretical values obtained by using a characteristic method to inte- 
grate an approximate form of the inviscid equations of motion (Ref. 4). 
Butler used the hypersonic slender body theory approximation to simplify the 
full equations of motions. Walkden and Caine estimated t ne pressure on the 
surface of a Butler wing at zero incident in a steady uniform stream by 
numerically integrating the two semi-characteristic forms of equations which 
govern the inviscid supersonic flow of an ideal gas with constant specific 
heat (Ref. 5). Squire obtained experimental results for a Butler wing at 
different Mach numbers and different angles of attach (Ref. 6,7). All 
previous numerical investigations have used inviscid equations and have 
excluded the wake region. In contrast, this study investigates the wake 

region as well as the effect of viscosity on the flow. The results of this 
study are compared with available experimental and numerical results. 

Graduate Research Assistant, Department of Mechanical Engineering and 
Mechanics, Old Dominion University, Norfolk, Virginia 23508. 

2 Eminent Professor, Department of Mechanical Engineering and Mechanics, Old 
Dominion University, Norfolk, Virginia 23508. 





In order to study the flow around a Butler wing, the Navier-Stokes 
equations are solved numerically. These equations are unsteady, compress- 
ible, viscous, and three-dimensional without neglecting any terms. The time 
dependency of the governing equations allows the solution to progress natur- 
ally from an arbitrary initial guess to an asymptotic steady state, if one 
exists. The equations are transformed from the physical coordinates to the 
computational coordinates, allowing the solution of the governing equations 
in a rectangular parallelepiped domain. The equations are solved by the 
McCormack time-split technique which is vectorized and programmed to run on 
the CDC VPS 32 (CYBER 205) computer. The codes were written in the 32-bit 
(half-word) FORTRAN which provides an approximate factor of two, decrease in 
computer time, and doubles the memory compared to the 64-bit word size. 
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2. GOVERNING EQUATIONS 

The governing equations for a thermal fluid system are the conservation 
of mass, momentum, and energy. These equations are developed for an arbi- 
trary region assuming the system is in continuum. Equations of motion for 
viscous, compressible, unsteady, heat conducting flow can be written as: 

Continuity: — + v • (p"u) =* 0, (2.1) 

at 

Momentum: + y . (puQ - T ) = o, (2.2) 


Energy: + v • (£u + <J-u*t)*0, (2.3) 

at 

V 2 

where E is the total energy per unit volume given by E = p (e + — + 

2 

potential energy + ...) and e is the internal energy per unit volume. 
Equations 2.3 can be simplified by assuming that the stress at a point is 
linearly dependent on the rate of strain (deformation) of the fluid (New- 
tonian fluid), 

3u, 3u. 8u. 

t.. = -P6..+u ( + — i ) + 6- ■ ] . (2.4 a) 

J J ax. ax. J ax, 

J i K 

where 6. , is the Kroriecker delta function, and y' is the second coeffi- 

\j 

cient of viscosity. The two coefficients are related to the coefficient of 
bulk viscosity (k) by the expression k = 2y/3+y‘. The contribution of k 
can be neglected if the pressure in a fluid is not changed abruptly during 
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its expansion or contractions. Under this assumption, the stress tensor can 
be related to the pressure and velocity components as: 


m i 




ij 


P 5 


8u i 3U . 


+ „ [( + — J - ) - - « i 

J 3x. 3x. 3 

J ' 


1J 


3u, 


3x, 


(2.4b) 
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This equation is valid for continuum flows. For an isotropic system, the 
heat flux in Eq. (2.3) can be expressed in terms of temperature gradient 
(Fourier's law of heat conduction) as: 
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q = - K 7 T 


(2.5) 


where K is the coefficient of thermal conductivity. 4 common approxima- 
tion used for viscosity is based on the kinetic theory of gases using an 
idealized intermolecul ar-forces potential. The relation is: 


(L .) 3/2 + $ o 


T 


T + S, 


(2.5) 


where 


S = 198.6* R 
o 

U r = 0.1716 np 


coefficients of thermal conductivity K can be determined from Prandtl 
number 



Pr 


(2.7) 
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where C y is the specific heat at constant volume and y is the ratio of 
specific heats. * 

It is necessary to have a supplementary relation to close the system of 
equations (Eq. (2.1) - (2.3)). By neglecting intermolecul ar forces (a ther- 
mally perfect systan), thermodynamic properties can be related as: 


P = pRT 


( 2 . 8 ) 


where R is the gas constant. Thermally perfect gas assumption allows 
expression of the internal energy (e) as a function of T only [e = e(T)]. 
in addition, assumption of calorically perfect gas [e(0) = OJ allows the 
following relation: 

e - C V T. (2.9) 

A combination of Eqs. (2.3) and (2.9) results in 

i 

i 

! P = pe (y~1) * (2J.0) 

I 

I . 

these equations (Eqs. (2.1) - (23)) are in conservative form. For sim- 
plicity, these equations can be expressed into a compact vector form as: 


a_u + a_F + _9 _g 
3t 3x 3y 


+ 



3Z 


( 2 . 11 ) 


where 
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puu - T xx + P 

PUV - T 


Eu + q x - 4» x + Pu 


puv - Tyx 

PVV - T + p 

y y 

PVW - T yz 

Ev + q y - + Pv 


PUW - T. 


PVW - T 


PWW > T + P 
Ew + q z - ij> z + Pw 


For the sake of generality, we can transform these equations from a pnysicai 
domain to a computational domain as 


at y 

5 , 


1- i _ 11 

35 35 35 


I y (2.12) 

J \ 8n an 8n / 


+ r 

•y 


3F aG an 


3c ac 3c/ 


The transformation coefficients can be computed from a functional re- 
lation between the computational coordinates and the physical coordinates. 


x = x(5,n,c), y = y(5,n,c), and z => z(5,n,c). 


(2.13) 
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5 = 5(x,y,z) 


(2.14) 


» n = n(x,y,z), and c = ?(x,y,z). 


If Eq. (2.14) is known, the transformation coefficients can be computed by 


direct differentiation. If the former relation is not known, after some 
algebraic manipulation, the transformation coefficients can be computed by: 
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3_z 

az 

35 

3n 

3C 



(2.15) 


where [J] and | J~ 1 1 are defined as (see Appendix A) 
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In the present case, the planes of grid are perpendicular to the x-direction 
thus allowing us to write: 


x = x(e) 

(2.17) 

y ~ y(n,c) 

(2. IB) 

z = z(n,c). 

(2.19) 


This reduces the metric coefficient from nine to five non-zero elements. 
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3. METHOD OF’ SOLUTION 


A time marching method is used to compute the solution. This allows us 
to capture the possible transient feature. This method is an explicit 
second-order accurate time-split predictor-corrector algorithm [3J. The 
governing equations (Eq. (2.12)) are discretized in computational direc- 
tions. In a compact form, they can be expressed as 


<j,k *[ L n‘“,,>] [Wl [ L C <4t t>] [ L t (1 V] [ L nK>l U ",j,k 0.1) 


where 


At = At 

n C 




and L._ , L , and L„ are the operators in e, n, and c directions, 

s n c 

respectively. A time step is completed in this algorithm with the appli- 
cation of each operator applied symmetrically about the middle operator. 
For example, operator can be defined as 


where 


L 5 (At ? ) = U 


out 

i 


(3.2a) 


Predictor step: 


U. . . = ui n . 

i,J,k i»J»k 


At. 

t. 

A? 


(F - F.^) — i + (G- - G ) M. i 
1 1 1 3x 1 3y 


8E 


+ (Hi - H m ) 

1 az 


J,k 


(3.2b) 


g 


5/ 








!«&>! 


Corrector step: 


1 (u? n . . +U. . 


M Olit 

U i ,j,k ‘ J r u i J,k 


At, 


AC' 


< F w- F i)?r 1 


(3.2c) 
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«,.i «W- "dr ' 


35 , 


O.k 
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This method has a time step stability limit, but there is no rigorous sta- 
bility analysis available for this. A conservative time step that is com- 
monly used is 


At < min 


Li!i + ill + + c 


Ax Ay Az 


JL + -L + -1 . 

Ax 2 Ay 2 Az 2 J 


1 “ 


(3.3) 


where c is the local speed of sound. 

In the supersonic region, there exists a large gradient which requires 
a very fine mesh to resolve it. If tney are not resolved, they produce a 
large oscillation which eventually blows up the ‘ elution. These oscil- 
lations of "low frequency" can be suppressed by adding a fourth order damp- 
ening. A coimion dampening used is the pressure dampening. ihis can be 
expressed in physical coordinates as 


-wi 


Kl + c 8 2 p 


36 . 


A? 


3 2 6 . 


3 U_ 

36 . 


I « 1,2,3 


(3.4) 


62 = n , and 13 
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where 6 1 = 5 , 



4. INITIAL AND BOUNDARY CONDITIONS 

In computational fluid dynamics the initial conditions usually corre- 
spond to a real initial situation for a transient problem, or a rough guess 
for a steady state problem. In practice, initial conditions are obtained 
from experiments, empirical relations, approximation theories, or previous 
computational results. An inappropriate initial guess may result in gener- 
ating unphysical ly strong transient waves which propagate through the compu- 
tational region dominating the flow field and eventually lead to a solution 
failure. In general, there are two important requirements that should be 
considered in the choice of initial conditions. First, they should be com- 
patible with the fixed upstream boundary conditions. Secondly, the initial 
conditions should be as physically close as possible to the actual nature of 
the flow field in the region under study. The former will minimize the 
lumber of iterations required for convergence. An attractive approach is to 
initialize the entire flow field (including the upstream boundary and the 
body surface) with a crude and simple guess (e.g., free stream condition). 
Then, during the course of the computation, both body and upstream boundary 
conditions are changed in a gradual manner to their final values over a 
prescribed number of iterations. The former approach is applied in only one 
step which is equivalent to impulsive initial conditions. 

It is equally important to implement a realistic, accurate, and stable 
method to determine boundary conditions. The application of certain condi- 
tions may cause numerical instability even though the flow is physically 
stable. There are neither mathematical nor physical justifications to im- 
plement a realistic boundary condition. Most of the boundary conditions 
currently implemented are drawn mainly upon intuition, wind tunnel experi- 
ence, and computational experimentation. There are three general types of 



boundary cond i tion:.. They are Dir.chlet conditions (specified function 
value), Neumann coiiditicus (specif i c -d normal gradient), and Robin condi- 
tions (a combination of both). Four important factors shouM be considered 
in tne election of boundary conditions. They are convergence, stability, 
computer time, and above all the physical justification. 

For this problem there are five different boundary conditions. They 
are upstream, downstream, lateral, top, and solid boundary. The upstream 
boundary conditions are the undisturbed free stream conditions and are lo- 
cated at a grid space away from the leading edge, i.e.. 


U 



upstream 


(4.1) 


A zero gradient in y-direction (parallel to the- primary direction of flow) 
is assumed for the downstream boundary, i.e., 


3 TT 
3y 


= 0 

downstream 


(4.2) 


The lateral boundaries are located far enough to avoid any influence on the 
interaction region. A boundary-! ayer profile can be prescribed on the lat- 
eral boundaries. These profiles can be obtained from their correspond ing 
points of a flow over a flat pate. Presently, a zero gradient in z-direc- 
tion is assumed for these boundaries, i.e., 

IE 

3z 


= 0 (4.3) 

1 ateral 


12 



The wall is assumed impermeable and no-slip boundary conditions are 
applied, therefore, all velocity components are assumed to be zero. The 
wall is also assumed to have a constant temperature Tw. A zero normal 
pressure gradient is assumed for the solid surface, i.e., 

if. 

3n 

This evaluation may appear to be based on the boundary-! ayer approximation 
(zero normal pressure gradient). In fact, it is a much milder approxima- 
tion, since constant pressure is not applied through the boundary layer but 
over one grid line in the boundary layer. This approximation has yielded 
stable computations for both non-separated and separated boundary layers 
[2]. From Eq. (A, 25b), Eq. (4.4) can be expressed as 


* 0 . (4.4) 

sol id 


3 P 
3n 


S n + £ n + c -n n + n 2 + n 2 + n 2 n 
x x y y z 7. P„ x y z P 
. _.v.. / — , ^ 1 n 

1/2 / . _ _ \ 1 / 2 


( n x + n y + n z) 

cn + c n + c ii 
^x x \v y ^z z 

/n 2 + n 2 + n 2 \ 

\ x y z ) 


n 2 + n 2 + n 2 


P r > 0 
1/2 5 . 


(4.5) 
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5. APPLICATION TO A BUTLER WING 

As mentioned in the Introduction, the Butler Wing is a delta wing which 
was proposed by 0. S. Butler [Ref. 4]. The planform of the body (Fig. la) 

is an isosceles triangle, and leading edges of the wing are laid along the 

Mach lines of the undisturbed stream. For the first 20% of the length, the 
body is a right circular core. The remainder of the body has elliptic sec- 
tions which become more eccentric as the sharp trailing edges are approached 

(Fig. lb). The semi major and minor axes are given by: 

Major axis (semi-span) = 

6 

)( 

Minor axis (thickness on = — 
centerline) 6 

x 

6 

where 3 2 = M 2 - 1 
00 

The wing is synmetric about (XZ) and (XY) planes. This permits us 
to use only one quarter of the entire physical domain (Fig. 1C) which is 
extremely advantageous from computational viewpoint. However, if the angle 
of attack is greater than zero then half of the physical domain should be 
considered. 

Some specific flowfield results nave been obtained for the Butler wing 
and these are discussed in the next section. 


0 < x < L (5.1a) 


0 < x < 0.2L (5.1b) 


1 - [ 


x - Q.2L, 
Q.SL ' 


0.2 < x < L 
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6. RESULTS AND DISCUSSION 


Grid generation is the very first-step which should be considered in 
obtaining flowfield solution over any configuration. Due to the data base 
management of our program, it is necessary to map the entire physical domain 
into a parallelepiped. Among the grid types, selection of an 0-type grid 
would produce a point singularity at the tip and a line singularity along 
the trailing edge. Nevertheless, this maps the solid boundary onto an 
entire face of the parallelepiped. A C-type grid is used in the present 
work. However, the solid boundary is not mapped onto an entire face of the 
computational box. This creates a potential problem in updating the bound- 
ary conditions. 

In this study, a two boundary grid generation (TBGG) technique develop- 
ed by R. E. Smith [1] is used. The method is essentially an algebraic 
method. The application of the TBGG method requires that the entire body be 
sliced into different cross-sections (Fig. 2a-b) . The cross-sections of 
wing are obtained in the stream-wise direction by analytical descriptions of 
the wing surface (Eq. 5.1). The TBGG method then is used to generate the 
grid for each station. There are fifty-five stations in the stream-wise 
direction, and each station has 64 x 36 grid points (Figs. 2C-2i). There is 
a total of 126,720 grid points which take 2.8 million 32-oit words of pri- 
mary memory (16 variables). The computational time required is 1.9 x 10~ 5 
sec/grid point/ iteration (2.5 sec/iteration). This is a typical require- 
ment for the CYBER 205 with two pipes. 

Results are obtained for a Butler wing at r*i w = 3.5, = 2 x 10 6 /ft, 
T w 3 390°R, Tw = 1092*R, L = 0.8 ft, and zero angle of attack. The com- 
puted pressures are plotted in Fig. 3. The pressure coefficient along the 
center line is shown in Fig. 3a. The results are compared with available 
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experimental and numerical results (Refs, 4,6,7). The results are in excel- 
lent agreement with the experimental results of Ref. 4. Moreover, they are 
closer to the experimental results than previous numerical results (Refs. 
4,6,7). Pressure ratios are plotted at constant chordwise position against 
the conical spanwise coordinates (41 .66% and 68.33X). These are seen to be 
in good agreement with experimental and numerical results (Fig. 3b-c). 
There are, however, same discrepancies in the results between 30° and 60°. 
This might be due to the fact that grids are not orthogonal near those 
regions. 

In the near future, extensive results will be obtained for different 
freestream conditions. These will be reported on in a future progress 
report. 
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7. CONCLUDING REMARKS 

General formulations are presented to investigate the flowfie'id over 
complex configurations for high-speed freestream conditions. An advanced 
algebraic method is used to generate grids around these configurations. The 
computational procedure developed is applied to investigate the flowfield 
over a Butler wing. Illustrative results obtained for specified freestream 
conditions compare very well with available experimental and numerical 
results. Further studies, however, are needed to establish the validity and 
versatility of the present code. After such model validations, it is anti- 
cipated to use the code tc investigate flowfield over complex configurations 
such as closed-bluff bodies (i.e., circular and elliptical cylinders on a 
flat plate, etc.). 
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Figure Zf, Grid distribution at Station Twenty-six 















Figure 2g. Grid distribution at Station Thirty-one 
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Figua 3a. Pressure 


coefficient along the center line 
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Fig. 3b. Pressure ratio at 41.67% of chord. 
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Figure 3c. 


Pressure ratio at 68.33% of chord. 



APPENDIX A 

MATHEMATICAL DETAILS FOR THE GOVERNING EOUATIONS 
A.l Curvil inear Coordinates 

In covariant coordinates system ( x ^ , the position vector of a point 
from the origin is expressed as 


r * Xj = ei xi + e 2 *2 * e3 } '3 


(A.l) 


where e. is the covariant base vector. 

In the present study, covariant coordinates are labeled x y, and z, 
and contravar iant coordinates are labeled n and c. The Cc iant base 
vectors are defined as 
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where J is the Jacobian of transformation. Magnitude of Jacobian (|J|) 
is the local value of the ratio of an elemental volume in the physical 
(usually distorted) cell to the corresponding elemental volume in the mapped 
(cubic) cell. 

The contravar i ant base vectors are defined as 
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Position vector can be expressed explicitly in terms of contravari ant 
vector (x*); however the infinitesimal vector dr can be expressed as 

dr a d x 1 = e i dx* ( A - 4 ) 

3x’ 

also the magnitude of arclength (ds) can be expressed as 

( ds) 2 adi" o dr»6 k ^ d x^ d x ^ («-5) 


where 6. is the kronecker delta, 

K 

6 . , = 1 if k = * 

kl 

6 k i =0 if k * *• 

Substitution of F.q. A. 4 into Eq. A. 5 will result in 
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dS2 = ( e i * e j)dx i dx. r g..dx i dx j . 

Wfl8re 9l 'J C3l1ed Covariant fundamental metric coefficients. 
These coefficients can be defined as 


9ij a C-H] ' [*H] .!£ 

3x, 


They are defined as 


911 3 x | + y| + z | 

912 * 921 ’ *5 \ * % * h *„ 

913 - 9 3 j • X £ X f ♦ y { y ? . ^ ^ 

922 3 * 2 + y 2 + z 2 

” n n 

923 a 932 a x n x c + y n y c + * n z c 
933 - x2 + y® + Z 2 


Similarly, covariant fundamental 


metric coefficients are defined 


g1J 3 V e i " CJ] [J] T =I!i 


9 n +<* +e® 

921 = 9 12 = ? n + r n + r 

x y y 2 n 2 

g!3 = q31 2 r 

9 'i Y C v + Cc+r- 

x X V s V ** 7 


n 2 + n 2 

y 2 


y y * v z 


(A.lUa) 




g 23 


(A.IOe) 
(A.IOf) ' 


“ 9 32 - * x C x + n y c y + n z c 2 
g 33 3 c 2 + c 2 + c 2 

Furthermore, there exists a unique relationship between contravariant 
and covariant fundamental metric coefficients, 


ij ^rs ^ls ” ^rt ^ls 


'J 


l9 (j l 


iSijl 


where 

C = g g . q 2 

11 22 33 23 


G = G =g g - g. g 

12 21 1323 12 33 


G 3 G =gg -g g 

13 31 12 23 13 22 


G = g g -g 2 

22 11 33 13 


G - G - g g-g g 

23 32 1213 23 33 


G = g g - g 2 

33 1 1 22 1.2 


(A. 11) 

(A. 12a) 
(A. 12b) 
(A. 12c) 
( A . 1 2 d ) 
(A.12e) 
(A.12f) 


There is also a relationship between covariant and contravariant base 
vector 



(A. 13) 
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wnere |g..j |J-1 1 3 

i .e. 


1 . ‘ Vt-Vn 1 1 - KWn) WVtV.I* 

|J-‘i 


e 1 = 


or 


e l 

|j-M 

(y n z cVn > 


( v { -Vn> 

e 2 


'w- Vs ! 

- |x sVVe> 

_e 3 _ 


_ (“eWe* 

- 1 \\-Vk ) 

(x eVV £ >_ 


(A. 14) 


where .] « e 


1 


There is also a relationship between contravariarit and covariant base 
vectors 


e j x e k 


li = 


|9 




where | g 1 | = | J { 2 


(A. 15) 


i.e. 


ei 

1 

"W 

< n yW : y> 

-fn C _r, C ) 
' x z Z X ' 

( n x c y*Vx) 


-(ty: z -e,c y ) 

(l r - r r ) 

^*’Z ’Z'x' 

■(s x c y *cy: x ) 

i3 


fyVW' 

-(5 n -K n ) 
' X z z *' 

(5 n -t n ) 
v z y y r 


or 


(A. 16a) 
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nr «C3»s 


[J- 1 ] « [ei e 2 e 3 ] 


<Vz-Vy> 

("xWx’ 

("x'/V** 


- I'/z-'z'y’ 
<'xV £ z c x> 

- (5 *'/- 5 y'x) 


(5 n -5 n ) 
' j z z y' 

-(£ n -5 n ) 
' X Z Z X’ 

(5 n -C n ) 
' z y y x 


(A. 16b) 


The relationship between vector bases can be obtained also by matrix alge- 
bra. From basic matrix identity, it can be written as 


p „ r , n Transpose of cofactor [J -1 ] . [[J* 1 ] ] 
[J] “ j^Tj " 


(A. 17) 


Equation (A. 17) is the same as Eq.(A.14). 


There exists an inverse relation for Eq. (A. 17), 


A. 2 Vector Representation in Curvilinear Coordinates 

A vector F(F^ i, F j, F z k) can be expressed in a contravarient co- 
ordinates as 


1/2 < 

F i - (g^) F 


(A. 18) 


where 


F 1 = lil F + ®*1 F y + F ? 
3 X x 3y y 3z 


(A. 19) 


Equation (A. 19) can be expanded as 




I 







A. 3 Normal Derivative in Curvilinear Coordinates 

A normal derivative of a scalar variable car be computed as 



where 


n 




VA = A i + A j + A K 

* jr Z 


or 


’ 

A i 


A 

X 

V 

' - [J] T 

5 

\ 

A 7 k 

1 

A 

z 


'1 


Equation (A. 23a) can be written as 


3A, 

„ [♦* *y ^ M 

fM 

3n 

(<> 2 f (j ,2 + ,j .2 )1 n 

A l 
n ( 


constant i|i y 1 

> V 


For constant - 5, Eq. (A. 24b) can be written as 


(A. 23a) 


(A. 23b) 
(A. 23c) 


(A. 24a) 


(A. 24b) 


(A.25«ah 


3A 

3n 


gll A + g 12 a 

3 t * n 3 c c 


(g 11 ) 


1 /2 


For constant - n. Eq. (A. 24b) can be written as 


3A 

3n 


g 21 A_ + q 22 A + g 23 A 
* 5 3 n 3 C 


For constant - c , Eq. (A. 24b) can be written as 


(A. 25b) 


also* 


3A 

3n 


g 31 A_ + g 32 A + g 33 A 
£ n C 




1/2 


(A); * 


le 1 ! 


V A 


3 

z 


Ip 1 ! ^ 


(e^*e J ) ki 


(A. 25c) 


/, 


e.g. 


JL J q ij L A . 


1/2 i= 1 


hu} J J ' 3 *' 


r I r 11 12 13 

(*>n “ — ,9 \ * 9 ft n + 9 ' V 


(A. 26) 


(A. 27) 
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A. 4 Miscellaneous Relations 

Angle between two grid lines is given by 


ft 


: > 
! 


4 2 


if 



1/2 


( A .28) ' 


/•/ 


7 ,..-r 


9, 


U 


COS 8 ,• j = — 

Ml 9jj! 


therefore for orthogonal grid, the following should be true, 


g.|j - 0 for i * j. 


(A. 29) 




Arclength is defined as 


fi 


\ '* 
f 


3 3 

(ds ) 2 = E r g . . dx 1 dx^, 
i=l j=l J 


(A. 30a) 


: V ; 


Arclength along x* coordinate is defined as 


- , i , .1/2 < 

(ds) * (g^) dx 


(A. 30b) 


The area of an element on which the x^ is constant is defined as 


(dr * | e d x"' x e d x* | = (g.,) dx^ dx^ 

■i "“u 1 J 


(A. 31) 


e.g. 




dr = | e 2 dn x e 3 d C. I = Gn drdd 


dr n - G 2 2 d€ d? 


(A. 32a) 
(A. 32b) 
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t <gr 



dr 3 £33 d £ d ri , 


Volume of an element is defined as 


dV |J-I| de dn dc 

d(C,n,c) 


1/2 

3 hi (e 2 xe 3 )j df, dn dc * ( |g i j | ) d? do dc (A. 33) 
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